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ABSTRACT 

N-Body simulations are an important tool in the study of formation of large scale structures. 
Much of the progress in understanding the physics of galaxy clustering and comparison with 
observations would not have been possible without N-Body simulations. Given the importance 
of this tool, it is essential to understand its limitations as ignoring these can easily lead to 
interesting but unreliable results. In this paper we study the limitations due to the finite size of 
the simulation volume. In an earlier work we proposed a formalism for estimating the effects 
of a finite box-size on physical quantities and applied it to estimate the effect on the amplitude 
of clustering, mass function. Here, we extend the same analysis and estimate the effect on 
skewness and kurtosis in the perturbative regime. We also test the analytical predictions from 
the earlier work as well as those presented in this paper. We find good agreement between the 
analytical models and simulations for the two point correlation function and skewness. We 
also discuss the effect of a finite box size on relative velocity statistics and find the effects 
for these quantities scale in a manner that retains the dependence on the averaged correlation 
function £. 

Key words: methods: N-Body simulations, numerical - gravitation - cosmology : theory, 
dark matter, large scale structure of the universe 



1 INTRODUCTION 

Large scale structures like galaxies and clusters of galaxies are be- 
lieved to have formed by gravitational amplification of small pertur- 
bations. For an overview and original references, see, e.g., Peebles 
(1980); Peacock (1999); Padmanabhan (2002); Bernardeau et al. 
(2002). Density perturbations are present at all scales that have been 
observed (Percival et al. 2007; Komatsu et al. 2008). Understanding 
the evolution of density perturbations for systems that have fluctu- 
ations at all scales is essential for the study of galaxy formation 
and large scale structures. The equations that describe the evolu- 
tion of density perturbations in an expanding universe have been 
known for a long time (Peebles 1974) and these are easy to solve 
when the amplitude of perturbations is much smaller than unity. 
These equations describe the evolution of density contrast defined 
as S(r, t) = (p(r, t) — p{t))/ p(t). Here p(r, t) is the density at r 
at time t, and p is the average density in the universe at that time. 
These are densities of non-relativistic matter, the component that 
clusters at all scales and is believed to drive the formation of large 
scale structures in the universe. Once the density contrast at rele- 
vant scales becomes large, i.e., \S\ > 1, the perturbation becomes 
non-linear and coupling with perturbations at other scales cannot be 
ignored. The equations that describe the evolution of density per- 
turbations cannot be solved for generic perturbations in this regime. 



N-Body simulations (Bertschinger 1998; Bagla and Padmanabhan 
1997b; Bagla 2005; Dolag et al. 2008) are often used to study the 
evolution in this regime. Alternative approaches can be used if one 
requires only a limited amount of information and in such a case ei- 
ther quasi-linear approximation schemes (Zel'dovich 1970; , 1989; 
Matarrese et al. 1992; Brainerd et al. 1993; Bagla & Padmanabhan 
1994; Sahni & Coles 1995; Hui & Bertschinger 1996; Bernardeau 
et al. 2002) or scaling relations (Davis & Peebles 1977; Hamilton 
et al. 1991; Jain et al. 1995; Kanekar 2000; Ma 1998; Nityananda & 
Padmanabhan 1994; Padmanabhan et al. 1996; Peacock & Dodds 
1994; Padmanabhan 1996; Peacock & Dodds 1996; Smith et al. 
2003) suffice. 

In cosmological N-Body simulations, we simulate a represen- 
tative region of the universe. This is a large but finite volume and 
periodic boundary conditions are often used. Almost always, the 
simulation volume is taken to be a cube. Effect of perturbations at 
scales smaller than the mass resolution of the simulation, and of 
perturbations at scales larger than the box is ignored. Indeed, even 
perturbations at scales comparable to the box are under sampled. 

It has been shown that perturbations at small scales do not 
influence collapse of perturbations at much larger scales in a sig- 
nificant manner (Peebles 1974, 1985; Little et al. 1991; Bagla & 
Padmanabhan 1997a; Couchman & Peebles 1998; Bagla & Prasad 
2008). This is certainly true if the scales of interest are in the 
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non-linear regime (Bagla & Padmanabhan 1997a; Bagla & Prasad 
2008). Therefore we may assume that ignoring perturbations at 
scales much smaller than the scales of interest does not affect re- 
sults of N-Body simulations. 

Perturbations at scales larger than the simulation volume can 
affect the results of N-Body simulations. Use of the periodic bound- 
ary conditions implies that the average density in the simulation 
box is same as the average density in the universe, in other words 
we ignore perturbations at the scale of the simulation volume (and 
at larger scales). Therefore the size of the simulation volume should 
be chosen so that the amplitude of fluctuations at the box scale (and 
at larger scales) is ignorable. If the amplitude of perturbations at 
larger scales is not ignorable then clearly the simulation is not a 
faithful representation of the model being studied. It is not obvious 
as to when fluctuations at larger scales can be considered ignorable, 
indeed the answer to this question depends on the physical quantity 
of interest, the model being studied and the specific length/mass 
scale of interest as well. 

The effect of a finite box size has been studied using N-Body 
simulations and the conclusions in this regard may be summarised 
as follows. 

• If the amplitude of density perturbations around the box scale 
is small (8 < 1) but not much smaller than unity, simulations un- 
derestimate the correlation function though the number density of 
small mass haloes does not change by much (Gelb & Bertschinger 
1994a,b). In other words, the formation of small haloes is not dis- 
turbed but their distribution is affected by non-inclusion of long 
wave modes. 

• In the same situation, the number density of the most massive 
haloes drops significantly (Gelb & Bertschinger 1994a,b; Bagla & 
Ray 2005). 

• Effects of a finite box size modify values of physical quantities 
like the correlation function even at scales much smaller than the 
simulation volume (Bagla & Ray 2005). 

• The void spectrum is also affected by finite size of the sim- 
ulation volume if perturbations at large scales are not ignorable 
(Kauffmann & Melott 1992). 

• It has been shown that properties of a given halo can change 
significantly as the contribution of perturbations at large scales is 
removed to the initial conditions but the distribution of most inter- 
nal properties remain unchanged (Power & Knebe 2006). 

• We presented a formalism for estimating the effects of a finite 
box size in Bagla & Prasad (2006). We used the formalism to es- 
timate the effects on the rms amplitude of fluctuations in density, 
as well as the two point correlation function. We used these to fur- 
ther estimate the effects on the mass function and the multiplicity 
function. 

• The formalism mentioned above was used to estimate changes 
in the formation and destruction rates of haloes (Prasad 2007). 

• It was pointed out that the second order perturbation theory 
and corrections arising due this can be used to estimate the effects 
due to a finite box size (Takahashi et al. 2008). This study focused 
specifically on the effects on baryon acoustic oscillations. 

• If the objects of interest are collapsed haloes that correspond 
to rare peaks, as in the study of the early phase of reionisation, we 
require a fairly large simulation volume to construct a represen- 
tative sample of the universe (Barkana & Loeb 2004; Reed et al. 
2008). 

In some cases, one may be able to devise a method to "cor- 
rect" for the effects of a finite box-size (Colombi et al. 1994), but 



such methods cannot be generalised to all statistical measures or 
physical quantities. 

Effects of a finite box size modify values of physical quanti- 
ties even at scales much smaller than the simulation volume (Bagla 
& Ray 2005; Bagla & Prasad 2006). A workaround for this prob- 
lem was suggested in the form of an ensemble of simulations to 
take the effect of convergence due to long wave modes into ac- 
count (Sirko 2005), the effects of shear due to long wave modes are 
ignored here. However it is not clear whether the approach where 
an ensemble of simulations is used has significant advantages over 
using a sufficiently large simulation volume. 

We review the basic formalism we proposed in (Bagla & 
Prasad 2006) in §2. We then extend the original formalism to the 
cases of non-linear amplitude of clustering and also for estimating 
changes in skewness and other reduced moments of counts in cells. 
This is done in §3. In §4 we confront our analytical models with 
N-Body simulations. We end with a discussion in §5. 



2 THE FORMALISM 

Initial conditions for N-Body simulations are often taken to be a re- 
alisation of a Gaussian random field with a given power spectrum, 
for details see, e.g., Bagla and Padmanabhan (1997b); Bertschinger 
(1998); Bagla (2005); Dolag et al. (2008). The power spectrum is 
sampled at discrete points in the k space between the scales cor- 
responding to the box size (fundamental mode) and the grid size 
(Nyquist frequency/mode). Here k is the wave vector. 

We illustrate our approach using rms fluctuations in mass 
a(r), but as shown below, the basic approach can be generalised 
to any other quantity in a straightforward manner. In general, a(r) 
may be defined as follows: 



a 2 (r) 



k 2tH y ' 



(1) 



Here P(k) is the power spectrum of density contrast, r is the co- 
moving length scale at which rms fluctuations are defined, k — 
\Jk 2 +ky + k 2 is the wave number and W(kr) is the Fourier 
transform of the window function used for sampling the density 
field. The window function may be a Gaussian or a step function 
in real or fc-space. We choose to work with a step function in real 
space where W{kr) = 9 (sin kr — kr cos kr) 2 / (fc 6 r 6 ), see e.g., 
§5.4 of Padmanabhan (1993) for further details. In an N-Body sim- 
ulation, the power spectrum is sampled only at specified points in 
the fc-space. In this case, we may write a 2 (r) as a sum over these 
points. 
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Here <r 2 (r) is the expected level of fluctuations in mass at scale r 
for the given power spectrum and a 2 (r, Lbox) is what we get in 
an N-Body simulation at early times. We have assumed that we can 
approximate the sum over the k modes sampled in initial conditions 
by an integral. Further, we make use of the fact that small scales 
do not influence large scales to ignore the error contributed by the 
upper limit of the integral. This approximation is valid as long as 
the scales of interest are more than a few grid lengths. 

In the approach outlined above, the value of a 2 at a given 
scale is expressed as a combination of the expected value o% and 
the correction due to the finite box size a 2 . Here a\ is indepen- 
dent of the box size and depends only on the power spectrum and 
the scale of interest. It is clear than cr 2 (r, L box ) < o"o( r )> an d, 
o\ (r, Lbox) > 0. It can also be shown that for hierarchical models, 
da 2 (r, Lbox)/dr < 0, i.e., af(r, Lbox) increases or saturates to a 
constant value as we approach small r. 

At large scales 0o(r) and a\ (r, Lb ox ) have a similar magni- 
tude and the rms fluctuations in the simulation become negligible 
compared to the expected values in the model. As we approach 
small r the correction term a\(r, L box ) is constant and for most 
models it becomes insignificant in comparison with 0o(r). In mod- 
els where <7 C 2 (r) increases very slowly at small scales or saturates 
to a constant value, the correction term a\ can be significant at all 
scales. 

This formalism can be used to estimate corrections for other 
estimators of clustering, for example the two point correlation. See 
Bagla & Prasad (2006) for details. 

The estimation of the rms amplitude of density perturbations 
allows us to use the theory of mass function and estimate a number 
of quantities of interest. For details, we again refer the reader to 
Bagla & Prasad (2006) but we list important points here. 



• The fraction of mass in collapsed haloes is under-estimated 
in N-Body simulations. This under-estimation is most severe near 
the scale of non-linearity, and falls off on either side. If we con- 
sider fractional under-estimation in the collapsed fraction then this 
increases monotonically from small scales to large scales. 

• The number density of collapsed haloes is under-estimated at 
scales larger than the scale of non-linearity. The maximum in col- 
lapsed fraction near the scale of non-linearity leads to a change 
of sign in the effect of a finite box-size for the number density of 
haloes at this scale: at smaller scales the number density of haloes is 
over-estimated in simulations. This can be understood on the basis 
of a paucity of mergers that otherwise would have led to formation 
of high mass haloes. 

• The above conclusions are generic and do not depend on the 
specific model for mass function. Indeed, expressions for both the 
Press-Schechter (Press & Schechter 1974) and the Sheth-Tormen 
(Sheth & Tormen 1999; Sheth, Mo, & Tormen 2001) mass func- 
tions are given in Bagla & Prasad (2006), and we have also checked 
the veracity of our claims for the Jenkins et al. (2001) mass func- 
tion. 



3 REDUCED MOMENTS 

In this section we outline how the formalism and results outlined 
above may be used to estimate the effect of a finite box-size on re- 
(2) duced moments. Reduced moments like the skewness and kurtosis 
can be computed using perturbation theory in the weakly non-linear 
regime (Bernardeau 1994). The expected values of the reduced mo- 
ments are related primarily to the slope of the initial or linearly 
extrapolated a (r), as all non-Gaussianities are generated through 
evolution of the Gaussian initial conditions and the initial o 2 (r) 
characterises this completely. We can use the expression for a 2 (r) 
as it is realised in simulations with a finite box size to compute the 
expected values of reduced moments in N-Body simulations in the 
weakly non-linear regime. 
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is the expected value of S3 for the given mode, i.e., when there 
are no box corrections and Ss 1 is the correction term in S3 due to a 
finite box size. Box size effects lead to a change in slope of a 2 , and 
hence the effective value of n changes. The last term is the offset 
in skewness in N-Body simulations as compared with the expected 
values in the model being simulated. We would like to emphasise 
that this expression is valid only in the weakly non-linear regime. 

In general we expect a\la% to increase as we go to larger 
scales. Thus the skewness is under estimated in N-Body simula- 
tions and the level of under estimation depends on the slope of 
o\ja\ as compared to the slope of a\. In the limit of small scales 
where a\ is almost independent of scale, we find that the correction 
is: 
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Here n is the index of the initial spectrum we are simulating. For 
non-power law models this will also be a function of scale. The 
correction becomes more significant at larger scales and the net 
effect, as noted above, is to under estimate S3. 

Similar expressions can be written down for kurtosis and other 
reduced moments using the approach outlined above. We give the 
expression for kurtosis below, but do not compute further moments 
as the same general principle can be used to compute these as well. 
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4 N-BODY SIMULATIONS 

In this section we compare the analytical estimates for finite box 
size effects for various quantities with N-Body simulations. Such a 
comparison is relevant in order to test the effectiveness of approxi- 
mations made in computing the effects of a finite box size. We have 
made the following approximations: 

• Effects of mode coupling between the scales that are taken 
into account in a simulation and the modes corresponding to scales 
larger than the simulation box are ignored. We believe that this 
should not be important unless the initial power spectrum has a 
sharp feature at scales comparable with the simulation size 1 . 

• Sampling of modes comparable to the box size is sparse, and 
the approximation of the sum over wave modes as an integral can be 
poor if the relative contribution of these scales to <7i is significant. 

Table 1 gives details of the N-Body simulations used in this 
paper. In order to simulate the effects of a finite box size, we used 
the method employed by Bagla & Ray (2005) where initial pertur- 
bations are set to zero for all modes with wave number smaller than 
a given cutoff k c . The initial conditions are exactly the same as the 
reference simulation in each series in all other respects. For a finite 
simulation box, there is a natural cutoff at the fundamental wave 
number kf = 27r/Lbox and simulations Al, Bl and CI impose no 
other cutoff. These are the reference simulations for the two series 
of simulations. Simulations A2, B2 and C2 sample perturbations at 
wave numbers larger than 2k f whereas simulations A3, B3 and C3 
are more restrictive with non-zero perturbations above Akf. The 
cutoff of 2k / and Akf corresponds to scales of 128 and 64 grid 
lengths, respectively. For the C series of simulations, the cutoff of 
2k f and Akf corresponds to scales of 80 h _1 Mpc and 40 h _1 Mpc, 
respectively. 

The background cosmology was taken to be Einstein-deSitter 
for the A and B series simulations. The best fit ACDM model from 
WMAP-5 (Dunkley et al. 2008) was used for the C series of simu- 
lations. 

In order to ensure that the initial conditions do not get a rare 
contribution from a large scale mode, we forced |5fej 2 = P(k) 
while keeping the phases random for modes k > 6k f . 

We have chosen to work with models where box size effects 
are likely to be significant, particularly with the larger cutoff in 
wave number. This has been done to test our analytical model in 
a severe situation, and also to further illustrate the difficulties in 
simulating models with large negative indices. 

We present results from N-Body simulations in the following 
section. 

4.1 Results 

We begin with a visual representation of the simulations. Figure 1 
shows a slice from simulations Al, A2 and A3 at two different 
epochs. The left panel is for the early epoch when r n i — 2 grid 

1 For example, simulations of baryon acoustic oscillations imprinted in the 
matter power spectrum may be affected by mode coupling even though the 
amplitude of fluctuations at the relevant scales is very small (Peebles 1974; 
Bagla & Padmanabhan 1997a; Takahashi et al. 2008). 



Model 


Description 






Cut Off(fc c ) 


Al 


Power Law, n — 


2.0 




k f 

J 


A2 


Power Law, n — 


2.0 




2k f 

J 


A3 


Power Law, n — 


2.0 




Akf 


Bl 


Power Law, n — 


2.5 




kf 


B2 


Power Law, n — 


2.5 




2k f 


B3 


Power Law, n — 


2.5 




Ak f 


CI 


ACDM, WMAP-5 BF, L box 


= 160 h" 


_1 Mpc 


kf 


C2 


ACDM, WMAP-5 BF, L box 


= 160 h" 


^Mpc 


2k f 


C3 


ACDM, WMAP-5 BF, L box 


= 160 h" 


^Mpc 


Ak f 



Table 1. This table lists characteristics of N-Body simulations used in our 
study. Here the spectral index gives the slope of the initial power spectrum 
and the cutoff refers to the wave number below which all perturbations are 
set to zero: kf = 2n/L hox is the fundamental wave mode for the simula- 
tion box. All models were simulated using the TreePM code (Bagla 2002; 
Bagla & Ray 2003; Khandai & Bagla 2008). 256 3 particles were used in 
each simulation, and the PM calculations were done on a 256' 1 grid. Power 
spectra for both the A and the B series of simulations were normalised to 
ensure it = 1 at the scale of 8 grid lengths at the final epoch if there is 
no box-size cutoff. A softening length of 0.1 grid lengths was used as the 
evolution of small scale features is not of interest in the present study. Sim- 
ulations for both the A and the B series were done with the Einstein-deSitter 
background and the C series used the WMAP-5 best fit (BF) model as the 
cosmological background, as also for the power spectrum and transfer func- 
tion. 



lengths in the model without a cutoff, and the right panel is for 
r n i — 8 grid lengths. The top row is for the simulation Al, the 
middle row is for the simulation A2 and the lowest row is for the 
simulation A3. The relevance of box size effects is apparent as the 
large scale structure in the three simulations is very different even 
at the early epoch when r n i = 2, much smaller than the effec- 
tive box size for these simulations. Disagreement between different 
simulations becomes even more severe as we go to the later epoch 
with r n i = 8 grid lengths. 

Visual appearance for simulations Bl, B2 and B3, shown in 
Figure 2 follows the same pattern. In this case the spectral index is 
closer to —3 than for simulations of the A series shown in Figure 1, 
hence the larger scale modes are more important for evolution of 
perturbations even at small scales. It is interesting to note that the 
largest under-dense region in simulation B 1 at early times is already 
comparable to the box size and hence we require L^ ox /r n i 3> 128 
for the effects of a finite box-size to be small enough to be ignored 
for simulations of the power law model with n — —2.5. This con- 
straint is even stronger for models with the slope of the power spec- 
trum closer to n = —3. 

Figure 3 shows the visual appearance for the C series of sim- 
ulations. Once again we find a significant change in appearance 
even with k c = 2k f at the earlier epoch, z = 1 in this case. 
This indicates that a box-size of 80 h _1 Mpc is insufficient if we 
wish to achieve convergence in the large scale distribution of mat- 
ter in models of cosmological interest. This reinforces conclusions 
of Bagla & Ray (2005) where we found that a box size of around 
150 h _1 Mpc is required for convergence in simulations of ACDM 
models. 
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Figure 1. The first, second and the third row in this figure show the slices for models Al, A2 and A3 (see table for details) respectively at an early epoch when 
the scale of nonlinearity is 2 grid lengths (left column) and a later epoch when the scale of nonlinearity is 8 grid lengths (right column) 
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Figure 2. The first, second and the third row in this figure show the slices for the models Bl, B2 and B3 respectively at an early epoch when the scale of 
nonlinearity is 2 grid lengths (left column) and a later epoch (right column) when the scale of nonlinearity is 8 grid lengths (right column) 
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Figure 3. The first, second and the third row in this figure show the slices for the ACDM simulations CI, C2 and C3 respectively at an early epoch (z 
(left column) and the present epoch (z = 0) (right column). 
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Figure 4. The first two rows in this figure show the average two point correlation function £ , and the next two rows show skewness S3 . The first and third row 
represent the early epoch and the second and fourth row represent the later epoch respectively. In all the panels models with k c = kf, k c = 2k j and k c = 4k f 
are represented by the solid, dashed and dot-dashed lines respectively. In all the panels, corresponding to every model in simulation (thick lines) theoretical 
estimates (thin lines) are also shown. Horizontal dashed lines in the lower rows shows the expected value of 53 in absence of any box-size corrections for the 
power law models. 
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Figure 5. This figure shows the mass function N(M)dM for the three series of simulations. The top row shows this for the early epoch and the lower row 
corresponds to the late epoch. In all the panels models with k c = kf,k c = 2k f and k c = 4k f are represented by the solid, dashed and dot-dashed lines 
respectively. In all the panels, corresponding to every model in simulation (thick lines) theoretical estimates (thin lines) are also shown. 



4.2 Clustering Amplitude 



an order of magnitude at some scales. Thus we may state that the 
model captures the essence of the box-size effects at large scales. 



The left column in Figure 4 shows the volume averaged correla- 
tion function f for the simulations being studied here. The top-left 
panel is for simulations Al, A2 & A3 at an early epoch (r ni = 2 
grid lengths in the model without a cutoff) and the second panel 
from top in this column shows f for the same simulations at a late 
epoch (r n i = 8 grid lengths. The corresponding plots in the second 
column show the same for the simulations in the B series, and the 
third column is for the C series of simulations. We have shown £ 
as a function of scale in these panels. Also shown are the linearly 
extrapolated values of f computed using our formalism for estimat- 
ing the effects of a finite box-size. Data from N-Body simulations 
is shown as thick curves whereas the theoretical estimate is shown 
as thin curves with the corresponding line style. It is clear that the 
analytical estimate for £ in a finite box captures the qualitative na- 
ture of the change from the expected values. The match is better at 
large scales where f is small and this is expected as the analytical 
estimate is linearly extrapolated whereas we are comparing it with 
results from an N-Body simulation. Our analysis works better for 
the n — — 2 model used in the A series of simulations and for the 
ACDM model in the C series of simulations as compared to the 
B series of simulations for the n = —2.5 model where it system- 
atically under-estimates the suppression of £. It is noteworthy that 
even in this case the differences between the simulation and the 
analytical model at large scales is of order of 20 — 30% whereas 
the box-size effect changes the clustering amplitude by more than 



4.3 Skewness 

The lower two rows in Figure 4 show the corresponding plots for 
5*3, shown here as a function of scale. Apart from the lines that 
show 53 from simulations (thick lines) and our analytical estimate 
for the weakly non-linear regime (thin lines), we also show the 
value of S3 expected in the weakly non-linear regime in absence 
of any finite box size effects for the three series of simulations. The 
analytical estimate of £3, computed using Eqn.(3) matches well 
with the values in N-Body simulation at large scales. It is notewor- 
thy that the match between the two is better for a larger cutoff in 
wave numbers. We believe that this is due to sparse sampling of 
the initial power spectrum at scales comparable to the box size and 
due to this our approximation of the sum over wave modes by an 
integral is not very good. 



4.4 Mass Function 

Figure 5 shows the number density of haloes for the three series of 
simulations as a function of mass of haloes. The haloes have been 
identified using the Friends of Friends (FOF) method with a linking 
length of 0.1 in units of the grid length. Plotted in the same panels 
are the expected values computed using the Press-Schechter mass 
function with a correction for the finite box size. 

For each series of simulations, and at each epoch, we fitted the 




Figure 6. The first two rows show the pair velocity as a function of distance, and the lowest two rows show the pair velocity as a function of average two point 
correlation function. The first and third row represent the early epoch and second and fourth row represent the later epoch. In all panels models with k c = kj 
(Al, Bl and CI), k c = 2kf (A2, B2 and C2) and k c = Akf (A3, B3 and C3) are represented by the solid, dashed and dot-dashed lines respectively. 




Figure 7. The first two rows show the pair velocity dispersion as a function of distance, and the lowest two rows show the pair velocity dispersion as a function 
of average two point correlation function. The first and third row represent the early epoch and second and fourth row represent the later epoch. In all panels 
models with k c = kf (Al, Bl and CI), k c = 2kf (A2, B2 and C2) and k c = Akf (A3, B3 and C3) are represented by the solid, dashed and dot-dashed lines 
respectively. 



12 Bagla, Prasad and Khandai 




Figure 8. The left panel shows the variation in the fractional correction in 53 i.e., Sz 1 /Sz (see Eqn. (4)), with the index of power spectrum at the scales 
Ltox/5 (top curve), L box /10 (middle curve) and L box /20 (lowest line). For a given tolerance of the error in S3 due to finite box effects, this gives us the 
largest scale at which the simulation may be expected to give reliable results. The middle panel shows contours of C\ at the scale of non-linearity (<ro = 1) 
for values Ci = 0.01 (top curve), 0.03 (middle curve) and 0.1 (lower curve). The contours are plotted on the L hoyL /r n i — (n + 3) plane and indicates the 
box size required for reliable simulations of a given model. The right panel shows contours of Sz 1 / Sz for the ACDM model that best fits the WMAP-5 data. 
Contours shown are for Sz^ / Sz = 0.01, 0.03 and 0.1. 



value of 8 C to match the simulation with the natural cutoff at the box 
scale. The same value of S c is then used for other simulations of the 
series. 

We find that the features of the mass function are reproduced 
correctly by the analytical approximation, namely: 

• The number density of the most massive haloes declines 
rapidly as the effective box size is reduced. 

• The number density of low mass haloes increases as the ef- 
fective box size is reduced. This feature is apparent only at the late 
epoch. 



4.5 Velocities 

In our discussion of analytical estimates of the effects of a finite box 
size on observable quantities, we have so far omitted any discussion 
of velocity statistics. The main reason for this is that the power 
spectrum for velocity is different as compared to the power spec- 
trum for density and one can get divergences for quantities anal- 
ogous to the second order estimators analogous to a 2 for models 
with — 3 < n < — 1. This is due to a more significant contribu- 
tion of long wave modes to the velocity field than is the case for 
density. Relative velocity statistics are more relevant on physical 
grounds and we use these for an empirical study of the effects of a 
finite box size on velocities. It is also important to check whether 
considerations related to velocity statistics put a stronger constraint 
on the box size required for simulations of a given model. 

We measure the radial pair velocity and also the pair velocity 
dispersion in the simulations used in this work. These quantities are 
defined as follows: 



(Kf) 

a 2 # 2 r 2 . 



(V) 



h(r) = - 



(6) 



where the averaging is done over all pairs of particles with separa- 
tion Tij = I (rj — Ti)\ = r. In practice this is done in a narrow bin 
in r. Here a is the scale factor, H is the Hubble parameter and Vi is 
the velocity of the ith particle. Similarly, the relative pair velocity 
dispersion is defined as: 



where v tj is the relative velocity for a pair of particles, and averag- 
ing is done over pairs with separation r. Dividing by a 2 _H r2 r 2 J gives 
us a dimensionless quantity and the usefulness of this is apparent 
from the following discussion. 

We have plotted the radial component of pair velocity as a 
function of scale r in the top two rows of Figure 6. Panels in these 
rows show the pair velocity for the different models at early and 
late epochs. In each panel, we find that the dependence of pair ve- 
locity on r is very sensitive to the small k cutoff used in generating 
the initial conditions for the simulation. It has been known for some 
time (Hamilton et al. 1991; Nityananda & Padmanabhan 1994) that 
h is an almost universal function of £. This is certainly true in the 
linear regime where h = 2£/3 for clustering in an Einstein-de Sit- 
ter universe (Peebles 1980). In order to exploit this aspect, and also 
to check whether the relation between h and £ in the weakly non- 
linear regime is sensitive to the box size, we plot h as a function of 
£ at the same scale in the last two rows of Figure 6. We find that all 
runs of a series fall along the same line and variations induced by 
the finite box size are small even in the non-linear regime. 

Figure 7 shows the relative velocity dispersion as a function of 
scale (top two rows) and also as a function of £ (lowest two rows). 
Again, we find that although the relative pair velocity dispersion at 
a given scale is sensitive to the size of the simulation box, its de- 
pendence on £ is not affected by the large scale cutoff. Thus we can 
use estimates of the correction in f to get an estimate of corrections 
in pair velocity statistics. 



5 SUMMARY 

The conclusions of this paper may be summarised as follows: 

• We have extended our formalism for estimating the effects of 
a finite box size beyond the second moment of the density field. 
We have given explicit expressions for estimating the skewness and 
kurtosis in the weakly non-linear regime when a model is simulated 
in a finite box size. 
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• We have tested the predictions of our formalism by comparing 
these with the values of physical quantities in N-Body simulations 
where the large scale modes are set to zero without changing the 
small scale modes. 

• We find that the formalism makes accurate predictions for the 
finite box size effects on the averaged two point correlation function 
£ and skewness. 

• We find that the formalism correctly predicts all the features 
of the mass function of a model simulated in a finite box size. 

• We studied the effects of a finite box size on relative velocities. 
We find that the effects on relative velocities mirror the effects on 
I 

It is desirable that in N-Body simulations the intended model 
is reproduced at all scales between the resolution of the simulation 
and a fairly large fraction of the simulation box. The outer scale up 
to which the model can be reproduced fixes the effective dynami- 
cal range of simulations. One would like S3 be within a stated toler- 
ance of the expected value at this scale. We plot Sz 1 /Sz for power 
law models at the scale Lbox/20, Lbox/10 and Lbox/5 in the left 
panel of Figure 8. These are plotted as a function of n + 3. It can 
be shown that this ratio, as also a\ / ao are functions of scale only 
through the ratio r/Lbox- We find that S , 3 1 /S , 3 is large for large 
negative n and decreases monotonically as n increases. This ratio 
is smaller than 10% only for n > 0.8 at r = Lbox/5. The corre- 
sponding number for r = Lbox/10 is n > —1.6, and for Lbox/20 
is n > —2.8 Clearly, the effective dynamic range decreases rapidly 
as n + 3 — > 0. This highlights the difficulties associated with 
simulating such models. 

Similarly, one would like a 2 and 0% to be comparable at the 
scale of non-linearity. From requirements of self similar evolution 
of power law models in simulations, we find that at the scale of 
non-linearity C\ < 0.03 is required for the effects of a finite box 
size to be ignorable. This gives us a lower bound on L hox /r n i for 
any given model. The middle panel shows the required Lbox/r n i 
as a function of n + 3 for Ci = 0.01, 0.03 and 0.1. Here Ci 
is the asymptotic value of a\ at r <C L box and is a fairly good 
approximation at small scales. We find that the required Lbox/r„; 
for n = —2 is more than 100 for C( = 0.03 at the scale of non- 
linearity. Thus we need a simulation with Lb ox > 10 3 if we are to 
probe the strongly non-linear regime (£ ^> 100) with some degree 
of confidence. Requirements for models with n < —2 are much 
more stringent, and for models like n = —2.5 even the largest 
simulations cannot be used to study the asymptotic regime. 

To put things in context for the favoured cosmological model, 
the right panel in Figure 8 shows contours of S^/S^ in the 
r - Lbox plane for the ACDM model that best fits the WMAP- 
5 data (Dunkley et al. 2008). We find that in order to ensure that 
the error in skewness is less than 10% at a scale of 10 h _1 Mpc, we 
need a simulation box of more than 200 h _1 Mpc. The required box 
size is much bigger if the tolerance on error in skewness is smaller. 
This is a very stringent requirement for simulations of the epoch of 
reionization where one would like to get the clustering right at the 
scales of a few Mpc. 

Given that the formalism we have proposed works well when 
compared with simulations, and the fact that calculations in this 
formalism are fairly straightforward, we would like to urge the cos- 
mological N-Body simulations community to make use of this for- 
malism. We would like to request simulators to report the fractional 
corrections to the linearly extrapolated amplitude of clustering and 
the fractional correction to skewness across the range of scales of 



interest. This will enable users of simulations to assess potential 
errors arising due to a finite simulation volume. 
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